Graphical Abstract

🚀 Please click on the images to enlarge them.

Reference-based compositional analysis

We used Kraken2/Bracken to perform taxonomic assignment of the raw sequences. The procedure included the remotion of reads mapping the rice genome (O. zativa, RefSeq accession number: GCF_000005425.2). The workflow is proposed by Lu et al. (2022), wrapped by Nextflow.

The analysis of the raw counts was developed by wrapping the output as a Phyloseq object to estimate effect of treatments and time over relative abundance, α-Diversity, β-Diversity and differentially abundant taxa.

MAG recovery

nf-core/mag was implemented to build and annotate the MAGs from the raw sequences. A co-assembly/co-binning strategy was chosen as an attempt to improve the quality and resolution of the MAGs. MEGAHIT was used for the assembly, while MaxBin2 and MetaBAT2 were the selected tools for binning. A de-replication and refinement step was carried by enabling DASTool within the workflow. For the downstream analysis, the curated MAGs were analyzed with MAGFlow/BIgMAG, that includes quality-measuring tools such as CheckM2, BUSCO, GUNC and QUAST, as well as the taxonomic annotator GTDB-Tk2 relying on database r220. The presence of Carbohydrate-Active enZYmes (CAZymes) was studied through dbCAN3, whilst Cluster of Orthologous Genes (COGs) were detected by COGclassifier. The selection of the MAGs was carried out using the criteria: completeness above 50% and contamination below 10%. In addition, MAGS belonging to either rice straw or soil with average differences <5% ANI were clustered using FastANI. For data visualization, Circos, iTOL and ggplot2 packages were used.

Fig. 1. Bioinformatics workflow followed to evaluate sequence quality, perform reference-based compositional analysis, and MAG recovery.

MAGs and CAZyme information

Fig. 2. a Two-dimensional Uniform Manifold Approximation and Projection (UMAP) embedding of normalized CAZy enzymes counts per MAG (before clustering with FastANI) colored by phylum (GTDB-Tk2). All CAZy classes (GH, GT, PL, CE, AA and CBMs), along with all their families, were considered for clustering using default values (15 neighbors and Euclidean distance). b Taxonomic classification (at the domain and phylum levels), quality (CheckM2), and identification as species or not of the clustered MAGs from the RS and soil metagenomes. The exact order of the missing phylum labels within the b plot can be read from the legend on a.

Fig. 3. Recovered MAGs from RS samples after clustering based on Average Nucleotide Identity (ANI). The visualization includes phylogenetic relationship among them (tree), specific highlighted phylum clades (colored sectors), genus classification (outer text), proportion of the main categories of Cluster of Orthologous Genes (COG) within each MAG (blue-orange-red-green pie chart), Carbon Active EnZyme (CAZyme) class normalized abundance (red scaled heatmap), classification as either high-quality or mid-quality MAG based on checkM2 measurements (blue binary), GUNC test passing (green binary), and genome size in Mbp (outer black bars).The displayed MAGs are representatives from the generated clusters using ANI, and they were selected based on the score: completeness – 0.5 x contamination. Conventions in CAZyme heatmap: GH glycoside hydrolases, GT glyco-syltransferases, PL polysaccharide lyases, CE carbohydrate esterases, CBMs carbohydrate-binding modules, AA auxiliary activities.

Fig. 4. Recovered MAGs from soil samples after clustering based on Average Nucleotide Identity (ANI). See Fig. 3 for feature explanation.

Fig. 5. Bubble plot depicting the number of normalized counts for CAZyme families associated with lignocellulose-degrading enzymes found in the MAGs recovered from RS samples. Visualization includes: MAG clustering based on normalized CAZyme counts of the selected families, and phylum information (GTDB-Tk2). The classification as MQ MAG or HQ MAG (CheckM2) is designated by the shape in the bubble plot. The displayed MAGs are representatives from the generated clusters using ANI, and they were selected based on the score: completeness – 0.5 x contamination. Displayed RS MAGs contain at least 30 normalized counts in 5 or more CAZyme families; complete information about normalized CAZyme counts in RS MAGs is found in Fig. S15.

Fig. 6. Bubble plot depicting the number of normalized counts for CAZyme families associated with lignocellulose-degrading enzymes found in the clustered MAGs recovered from soil samples. See Fig. 5 for feature explanation. All soil MAGs are displayed.

Fig. 7. Boxplot showing the distribution and abundance of normalized counts for CAZyme families associated with lignocellulose-degrading enzymes found in the MAGs recovered from RS and soil samples.

β‑Glucosidase activity

Fig. 8. Left: bubble plot depicting normalized CAZyme counts for β-glucosidase families (GH1, GH3, GH5 and GH9) per sample; GH glycoside hydrolases. Right: data distribution of experimental β-glucosidase activity per sample. The displayed normalized CAZyme values represent the maximum count in a given MAG of each GH family within the sample; treatment colors are shared by both the bubble plot and the boxplot.

Field experiments

Fig. S1. Field diagram of plots and yard distribution, illustrating yard sizing and color plotting distribution on the field. Green dot lines indicate rice furrows without treatment used as a natural barrier between plot treatments. Blank spots represent spaces without treatments due to field irregular terrain impeding a treatment application. Arrow indicates magnetic north.

Fig. S2. Illustration of sampling points for soil and straw for physicochemical and enzymatic analysis on field. A shows the distribution of soil sampling points within each plot; B depicts the distribution of rice straw sampling points within each plot; C displays a field picture showing the marks used for the replicas and color code per treatment: pink represents Cont plots, blue FBa plots, yellow FBaNit plots, and violet FCent plots; D presents RS manual collection; E exemplifies a tube probe for soil sampling; and F,G exhibit soil sample disposal in marked Ziploc® bags.

Soil parameters

Fig. S3. Experimental measures of physicochemical parameters in soil samples including total carbon, total nitrogen, pH and phosphorus. The data were normalized using the Min-Max method to adjust the raw values to the same scale; the displayed numbers represent the mean of the quantified parameter per sample.

Enzymatic activity

Fig. S4. Experimental protease, acid, and alkaline phosphatase activities of RS (a) and soil (b) samples. The enzymatic activity was measured according to the guidelines established by Alef & Nannipieri (1995).

Reference-based analysis

Fig. S5. Screenshot of a Krona plot displaying the taxonomic classification and relative abundance, performed with Kraken2 coupled to Bracken, of the reads belonging to the mock community sequenced as part of quality control for the DNA extraction process.

Fig. S6. Taxonomic profiling of the sequences using Kraken2 coupled to Bracken: a Relative domain abundance per matrix source, merging all the replicates per sample. b α-Diversity, colored by treatment, measured through Shannon and Inverse Simpson indices. c Sample ordination through Principal Coordinate Analysis (PCoA) using a Bray-Curtis dissimilarity matrix as input. The ellipses represent 95% confidence regions around the centroids of each group calculated by the assumption of multivariate t-distribution. d Two-dimensional Uniform Manifold Approximation and Projection (UMAP) embedding of β-Diversity estimation of the samples based on the Aitchison distance. e Compositional variability of the samples, measured as distance to the centroid in a multivariate space computed with the Aitchison distance among samples. ANOVA analysis results show no significant differences in terms of group dispersion. f Proportion of the variance explained by the factors Time and Treatment, segmented by matrix source, in Permutational Analysis of Variance (PERMANOVA) implemented according to different distance/dissimilarity metrics; the significance of each factor is determined by the symbol color. g Selected differentially abundant genera in RS samples based on Analysis of Compositions of Microbiomes with Bias Correction 2 (ANCOM-BC2). Left: log₂FC of treatments to Cont0, and Right: Boxplots of selected genera abundance (counts) in log₁₀ scale. All displayed genera exhibit p-value < 0.05 adjusted with Holm–Bonferroni method for multiple comparisons. Cont0 and Cont1 in RS comparisons indicate Cont samples at time 0 and 1, respectively. Differences throughout time for all treatments against Cont samples at time 0 are shown. All the analysis from b to f were carried out only with Bacteria data and normalized counts as Counts per million (CPM). The treatment colors are transversal to plots b, c, d, e, f and g. ANCOM-BC2 was implemented with agglomerated data at genus level, and with the most abundant taxa (> 0.5%).

Fig. S7. Percentage of annotated reads using Kraken2/Bracken segmented by matrix and time point. For RS at time 0 samples were taken from a bulk of material; this same bulk was subsampled with specific amounts of rice straw to start the experiments, and hence accounting with only 1 sample at time 0.

Fig. S8. Relative phylum (a, c) and genus (b, d) abundance of agglomerated data (Bacteria, top; Eukaryota, bottom). Displaying abundance values that belong to the portion of the reads annotated in each domain (see Fig. S6a and Fig. S7). Replicates were merged, and the panels are split by matrix source. On phylum plots, only phyla with relative abundance above 0.5% are displayed, whilst on genus plots genera above 3% and 2% are depicted for Bacteria and Eukaryota, respectively.

Fig. S9. Relative genus abundance of agglomerated RS data with panels split by treatment. Only genera with relative abundance above 3% are displayed.

MAG reconstruction

Fig. S10. Metrics summaries of the MAGs built from RS (a) and soil (b) samples; the columns represent, in order, the proportion of: mid-quality and high-quality MAGs according to CheckM2 estimations; MAGs passing the GUNC test; and annotated MAGs at different taxonomic levels based on GTDB-Tk2.

Fig. S11. Presence/absence of MAGs annotated with GTDB-Tk2 in RS samples at genus level. Genome size, in Mbp, is indicated as the point size, while the color scale represents the relative abundance of each MAG per sample, and the point shape shows whether the displayed MAG is classified as high-quality or mid-quality (CheckM2). When MAGs with the same genus annotation were found in the same sample, the most complete one is depicted and their abundances are added up.

Fig. S12. Presence/absence of MAGs annotated with GTDB-Tk2 in soil samples at genus level. Genome size, in Mbp, is indicated as the point size, while the color scale represents the relative abundance of each MAG per sample, and the point shape shows whether the displayed MAG is classified as high-quality or mid-quality (CheckM2). When MAGs with the same genus annotation were found in the same sample, the most complete one is depicted and their abundances are added up.

MAGs and CAZyme information

Fig. S13. Circos plots showing the proportion of the 4 main categories (Metabolism; Information Storage and Processing; Cellular Processes and Signaling; and Protein Coding Sequences with Unknown function) of Clusters of Orthologous Genes (COGs) found within the MAGs recovered from RS (top) and soil (bottom) the samples. The inner ring accounts for the total number of COGs related to every sample or COG category; the outer ring depicts the relative abundance of COGs from each sample or COG category; the width of the ribbons linking any sample and a COG category depicts their relative abundance to each other. The number of counts per category was normalized using the total number of genes found per MAG.

Fig. S14. Two-dimensional Uniform Manifold Approximation and Projection (UMAP) embedding of normalized CAZy enzymes counts per MAG (before clustering with FastANI) colored by treatment and shape by the time in which the MAGs were recovered.

Fig. S15. Clustered MAGs recovered from RS samples and their normalized counts for CAZy families associated with lignocellulose-degrading enzymes.The organization of the tree relies on the phylogenetic information retrieved from GTDB-Tk2. The displayed MAGs are representatives from the generated clusters using ANI, and they were selected based on the following score: completeness – 0.5 x contamination.

GTDB taxonomy of RS MAGs

GTDB taxonomy of soil MAGs